Textures and non-Abelian vortices in atomic d-wave paired Fermi condensates 



o 
o 

(N 



H. M. Adachi 1 , Y. Tsutsumi 1 , J. A. M. Huhtamaki 1 * 2 , and K. Machida 1 
1 Department of Physics, Okayama University, Okayarna 700-8530, Japan and 
2 Department of Applied Physics, Helsinki University of Technology, P.O. Box 5100, 02015 TKK, Finland 

We report on fundamental properties of superfluids with d-wave pairing symmetry. We consider 
neutral atomic Fermi gases in a harmonic trap, the pairing being produced by a Feshbach resonance 
via a d-wave interaction channel. A Ginzburg-Landau (GL) functional is constructed which is 
symmetry constrained for five component order parameters (OP). We find OP textures in the cyclic 
phase and stability conditions for a non-Abelian fractional 1/3- vortex under rotation. It is proposed 
how to create the intriguing 1/3- vortex experimentally in atomic gases via optical means. 

PACS numbers: 03.75.Ss, 67.85.-d, 03.75.Mn 



S3 ' 
O 1 

o , 
i 

s : 

a ■ 



-a 
c 

o 
o 



> 

(N 
l> 
On 
(N 

l> 
O 

&\ 
O 



X 



Superfluids with s-wave pairing symmetry have been 
realized by using a Feshbach resonance of 6 Li atom gases 
at i?=822G in 2OO50. It is natural to expect that the 
research front of cold atom gases develops towards realiz- 
ing condensates with higher relative angular momentum 
of the Cooper pairs. In fact, much effort from theoretical 
and experimental sides is now focused on p-wave pairing 
in 6 Li where ap -wave Feshbach resonance at H—159G is 
already confirmed to exist Q • Although p-wave molecules 
have been formed 0, experimental evidence of p- wave su- 
perfluidity in 6 Li has not reported. 

At this stage it might be useful to theoretically investi- 
gate d-wave superfluidity in neutral Fermi gases in order 
to further motivate experimental and theoretical works 
towards this direction. Hulet has performed a coupled 
channels calculation and found a Feshbach resonance for 
d-wave channel in 6 Lii, and hence we have a good rea- 
son to explore this possibility from a theoretical point of 
view. 

It is known that the pairing symmetry of high T c 
cuprates is described by the d x i_ y i state. The study 
of superfluids with d-wave symmetry is important be- 
cause several strongly correlated superconductors, or so- 
called heavy Fermion superconductors belong to uncon- 
ventional pairing states, such as d, d + id, or /, etc, in- 
cluding high T c whose pairing mechanism is still unclear. 
It is also interesting because of the richer physics asso- 
ciated with the many internal degrees of freedom of the 
relative orbital angular momentum I > 1 of the Cooper 
pairs compared to the I = (s) and I = 1 (p) cases. In 
particular, the spatial structure of the order parameter 
(OP), or the textures and vortices, which are hallmarks 
of superfluidity, are intriguing in the confined ultracold 
atomic gases. 

We note here that the OP structures of the <i-wave 
superfluid of fermionic atoms and spinor condensates of 
bosonic spin-2 atoms are mathematically very similar be- 
cause both OP's have 5 components. Whereas the lat- 
ter has been extensively investigated already [3, HI, the 
former has not been studied so far in connection with 
neutral atom gases. As we will see soon, the boundary 
conditions due to the harmonic trap are different between 



the two cases in an essential way, giving rise to novel tex- 
tures and vortices in the former. This difference arises 
because of the different origin of the OP degeneracy due 
to the internal degrees of freedom: The orbital angular 
momentum lives in real space in the former unlike the 
spin degrees of freedom in the latter although both are 
of the same 50(3) symmetry in a homogeneous system. 

Previously, in his seminal work, Mcrniin (i presents 
a general framework based on Ginzburg-Landau (GL) 
formalism to describe a d-wave superfluid in an infinite 
bulk system. He exhaustively classifies the ground state 
phase diagram into three regions; ferromagnetic (FM), 
polar (PO) and cyclic (CY) phases, depending on the 
coupling constants or fourth order coefficients in the GL 
functional. The cyclic phase, which we focus on in this 
paper, is the most intriguing one because fermion super- 
fluids with s-wave (I — 0) or spinless p-wave (/ = 1) pair- 
ings produced by a Feshbach resonance via a magnetic 
sweep, or the spin-1 spinor BEC0, HI do not support 
this phase which only appears above I > 2. 

The main motivations of this paper are to provide the 
fundamental physical properties of a d-wave superfluid 
confined by a harmonic potential in two-dimensions (2D) 
in order to help to identify the d-wave nature experi- 
mentally and to present a non-Abelian 1 /3- vortex which 
is the energetically favored state under external rota- 
tion. The existence of a 1/3- vortex provides a spectro- 
scopic means to characterize d-wave superfluidity. Mu- 
tually non-commutative 1/3- vortices themselves are al- 
ready pointed out to be allowed topologically in an F = 2 
spinor BEC by several authors[lO(. However, there is no 
serious calculation to examine its stability from energetic 
point of view neither in spinor BEC nor in d-wave fermion 
superfluids, which is one of our main purposes in this pa- 
per. 

The non-Abelian 1 /3- vortices might turn out to be use- 
ful for quantum computation because a spatial arrange- 
ment of mutually non-commutative 1/3- vortices can be 
used for storing information, namely they can be used as 
a topologically protected qubit. This is somewhat similar 
to the idea which utilizes the Majorana particles bound 
to the vortex core in chiral p + ip superconductors 11 1 . 
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Note that the 1 / 3- vortex discussed here does not support 
Majorana zero energy particles. 

Our arguments are based on GL theory which is gen- 
eral enough to describe a d-wave superfluid in terms of 
the free energy expanded up to fourth order in the OP. 
GL theory is valid near the transition temperature T c , 
however, its range of applicability is known to be wider 
empirically. The OP A for a d-wave superfluid is spanned 
by the spherical harmonics Y\ m (k) (k is a unit vector on 
the Fermi sphere) 



A(r) = A m (r)Y h: 



(k) = kiBykj 



(1) 



with I = 2, m = -2,-.- ,2 and i,j = 1,2,3. The re- 
peated indices are summed over. The coefficients A m (r) 
are a complex valued functions of r = (a;, y). To empha- 
size the essential points, we consider a two-dimensional 
system assuming that the obtained objects extend uni- 
formly to the third dimension. Alternatively we may use 
the symmetric traceless 3x3 matrix Bij. The i?-matrix 
which has five independent elements is convenient in con- 
structing the GL functional for the OP with 50(3) sym- 
metry in addition to U(l) gauge symmetry. In the follow- 
ing we use these two notations interchangeably. The bulk 
GL free energy functional fbulk is derived by Mermin as 



The other cyclic states are obtained from this by rota- 
tions. The weak coupling estimate for the Fermi sphere 
(}3 2 = 2 fa and /3 3 = 0) predicts that the stable phase is 
on the boundary between the CY and FM phases. 

The gradient terms are constructed by considering the 
possible contractions of (diBj k )* and (diB mn ); namely 
there are three independent terms: (1) (diBj k )* (diBj k ) 
(2) (diB^YidkB^) and (3) (c^)*^^). This can 
be understood from the fact that the decomposition of 
the two angular momenta (L = 1) X (L = 2) — ► L = 3, 2, 1 
where the former (latter) corresponds to d (OP). Then, 
the gradient energy fill, fl^ can be summed up as 



farad = Ki(diB jk )*(diB jk ) + K^diB^idkBij), 



K^djBuYtdkBik) 



(4) 



J bulk 



aTrB*B + (3i\TrB 
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where Ki = K, K 2 = 2K, K 3 = 2K with K = 2if /35 in 
the weak coupling approximation which also gives a — 
2a /15, /3 2 = 8/3 /315, OP amplitude A = ^JWv 
and the coherence length £o = \J Kq I &o ■ (For o-o, /3oj 
and K , see [bj]). We can derive the result of Eq. (4) 
also starting from the A m representation after simple but 
tedious calculations. 

In addition to these bulk and gradient terms, we take 
into account a harmonic trap potential: 



/3 2 (TrB*By + (3 3 Tr(B 2 B* 2 ) 



(2) 



where a(T) — aa(T — T c ) and T c is the transition temper- 
ature. There are three independent fourth order terms 
/?2 and /?3. The trace operation for a matrix is de- 
noted by Tr. This bulk energy fbulk is recast into 



fbulk = a (T - TcjY.^Ai 
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External rotation may be included by replacing d — > d — 
i^Pn x r where the rotation axis fl || z. Therefore the 
resulting total GL functional is given by 



f total — fbulk fgrad f harmonic 4~ fcent 



(6) 



where fcent is the centrifugal potential [3]. The associ- 
ated mass current is derived as 



+ (/3i + g/3 3 )|e| 2 --/3 3 |/T}S l |^| 4 (3) 

where the orbital singlet pairing amplitude O = 
(—l)' l AiA-i/\Ai\ 2 and the orbital momentum / = 
A*FijAj/\Ai\ 2 with F being the 5 x5 spin matrix The 
three phases are characterized by FM ((©) = 0, (/) ^ 0), 
CY ((6) = 0,</) = 0), PO «0) 0,(/) = 0) 
where (• • • ) is the ground state average. It is clear 
from Eq. (3) that CY is stable when /3i + \fc > 
and /?3 < 0, while FM and PO phases occupy the 
other regions of the (03, 0i) parameter space (see Fig. 
1 in Ref. Q). The canonical CY phase is of the form 
A(k) = iF 22 (k) + V2F 20 (k) + iY2-2(k) or in a vecto- 
rial notation: (A 2 , Ai, Aq, A-i,A- 2 ) T = (i, 0, y/2, 0, i) T . 



ji = 21m[KiB* k ViB jk + K 2 B* k V k B ij + K 3 B^kB jk ]. 

The following numerical calculations arc done in two- 
dimensional plane of cartesian coordinates (x,y). The 
coupled GL equations for five components are solved on 
discretized lattice using 101x101 meshes via a simple it- 
eration. The OP amplitudes and length are scaled by Ao 
and £o respectively. The trap frequency u> is normalized 
by 1 h J 2 J? with ui the atomic mass. We fix u> = 1/25 in 
these units throughout this paper. 

Before discussing the CY and 1/3- vortex, we briefly 
touch on the weak coupling case, which indicates the sta- 
ble phase on the boundary between FM and CY for an 
infinite system as mentioned, in order to understand the 
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FIG. 1: (Color online) Stable ferromagnetic texture (FM) in 
the weak coupling case at rest. The f x -fy vectorial pattern is 
shown. The length is scaled by £o- 



importance of the boundary condition through the gra- 
dient terms. We have done numerical calculations for 
a trapped system by solving the GL equations for five 
components. It turns out that FM is energetically favor- 
able over CY. This is because of the boundary condition, 
which strongly constraints the possible phase. In fact, 
FM is more flexible than CY in the sense that the /- 
vectors can arrange to stabilize the state. As seen from 
Fig. 1 where the /-vector configuration is shown for FM, 
the /-vectors tend to align along the circumference. This 
arrangement is energetically advantageous because the 
four point nodes situated along the / direction tend to 
lie outside the system, leading to maximal gain in the 
condensation energy. This is analogous to that in p-wave 
superfluids|l4j. This flexible feature is absent in CY be- 
cause / = 0. Thus under confinement the degeneracy 
between FM and CY is lifted, the latter being never sta- 
bilized over FM in the weak-coupling limit. 

Let us now come to our main part. We take the Pi val- 
ues appropriate for CY, for example, P1/P2 = L P3/P2 — 
— 1. The following results do not depend on these val- 
ues and are generic in the cyclic region. By solving 
the GL equations, we obtain the stable CY in the pres- 
ence of a trap. Among various CY forms derived from 
the canonical CY we stabilize the particular CY phase, 
called CY-z described by (1, 0, 0, a/2, 0) T , which is ob- 
tained by exp(i cos~ 1 (^)F y ) exp(-z|F 2 )(i, 0, y/2, 0, i) T . 
This CY-z will be seen to support the 1/3- vortex later. 
As seen from the main panel of Fig. 2 where the cross sec- 
tions of each component are displayed, CY-z is dominant 
in the central region because (f ) = and simultaneously 
(0) = as shown in Figs. 2(a) and (b). The surface 
region consists of FM or PO phases, which are intermin- 
gled. Note that FM is advantageous in the surrounding 
region because as mentioned the system can gain in the 
condensation energy by tuning the direction of /. This 
OP texture differs completely from the spinor F=2 BEC 
where the CY-z extends all the way to the surface of the 
cloud[l5| because of different boundary conditions. We 
emphasize that wherever OP varies spatially, the gradient 




FIG. 2: (Color online) Stable texture of the cyclic phase (CY- 
z) at rest. Cross section is shown along the horizontal direc- 
tion together with their contour maps of the five components 
as seen from above. |/| (a) and |0| (b) are displayed. The 
field of view is 12£o x 12£o. 

coupling terms f gra d m Eq. (4) inevitably induce other 
components. 

The gap structure of CY consists of 8 point nodes; In 
canonical CY (i, 0, \/2, 0, i) T the point nodes are situated 
at the 8 directions (±1, ±1, ±1) given by 8 corners of the 
inscribed cube inside the Fermi sphere. In CY-z the in- 
scribed cube is rotated so that the (1,-1,1) direction lies 
now parallel to the z axis. Thus this CY-z is three-fold 
symmetric with respect to combined gauge transforma- 
tion and rotations about the z axis. This is the origin 
for the possible existence of a 1/ 3- vortex in CY-z as dis- 
cussed below. The mass current given above flows spon- 
taneously circularly around the center of the trap and 
gradually decreases farther away from the origin. 

By rotating CY-z we can create the 1/3-vortex whose 
main structure is described by {e l6 , 0, 0, V2, 0) T in the 
central region as seen from Fig. 3 where the cross section 
for each component is displayed. Note that at the center 
the m = +2 component vanishes whereas the m = — 1 
component is non-zero. Simultaneously the other com- 
ponents are induced around the surface region. Com- 
bination of the three-fold symmetry around the z-axis 
for CY-z with an additional gauge transformation leads 
to the 1/3-vortex form. Namely, CY-z Y 2 ,2 + V^Y2,-i 
transforms into Y^e 2 ^ + v / 2Y2,_ie -4 * under a rotation 
of an angle (f> around the z-axis, which can be rewritten 
as e-^(r 2 ,2e 3 ^ + V2Y 2 ,-i) = e- l9 / 3 (Y 2 , 2 e ie + \/2Y 2 _i) 
after identifying cf> — #/3. This is apparently of the 
1/3-vortex form (e l6 , 0, 0, 0) T . The vortex center is 
dominated by the j4_i component where the A 2 compo- 
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FIG. 3: (Color online) The 1/3- vortex is stable at fi = 0.20w. 
Main panel shows the cross section for five components along 
the horizontal direction together with their contour maps as 
seen from above. |/| and \Q\ are displayed in (a) and (b), 
respectively. The field of view is 12£o X 12£o- 
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FIG. 4: (Color online) Stable texture at high rotation at 
Q, — 0.28ij;. I/I (a) and |6>| (b) are displayed where the cen- 
tral region is a 1/3- vortex surrounded by different kinds of 
vortices. The field of view is 20£o x 20£o- 



It is remarkable that this 1/3- vortex is the most stable 
vortex among various possible ones which we have solved 
in order to compare their energies. The critical rota- 
tion frequecy is found to be Q cr — 0.24u> at T = 0.5T C 
beyond which a single 1/3- vortex becomes energetically 
stable compared to the CY-z. Upon increasing fl fur- 
ther, different kinds of vortices emerge from the surface 
as shown in Fig. 4. They are similar to the half-quantum 
vortex (HQV) seen in p-wave supernuids[T4j. At the cen- 
ter the 1 /3- vortex exists surrounded by 8 HQV like ob- 
jects. In order to create non-commutative 1/3-vortices in 
a system and to observe non-Abelian braiding statistics 
among them, we need a special devise, for example, after 
creating a 1/3- vortex we should change the rotation axis 
from the z-axis to a different one. 

The initial preparation of CY-z itself is not difficult 
where the appropriate atom population ratio (1:2) in A 2 
and A— 1 is needed. Since these two components are dif- 
ferent orbital angular momentum states, the population 
in the components can be adjusted using Raman tran- 
sitions, similarly as in F=2 spinor BEC[16[ or by using 
RF transition in the presence of Zeeman splitting of the 
5 components under an applied field as was done by Hi- 
rano group [171]. Furthermore, once CY-z is prepared, it is 
possible to imprint the unit phase winding in the A 2 com- 
ponent using a circular polarized Raman transition (l6j. 
It may be difficult to use the optical spoon method to 
create the vortex, which was utilized to produce vortices 
in scalar BECsflg. 

We thank T. Ohmi, M. Ozaki, T. Mizushima, M. 
Kobayashi, T. Kawakami, T. Hirano and R. Hulet for 
useful discussions. 

During preparing the manuscript, we became aware of 
related preprints [ij], [2(|. The former (latter) treats d 
(/)-wave pairing in atomic gases from different points of 
view than here. 



nent vanishes due to the phase singularity at the core 
as seen from Fig. 3. Thus the core of the 1/3- vortex 
is ferromagnetic. The nodal structure of the ferromag- 
netic core region is characterized by a line node on the 
equator of the Fermi sphere in addition to point nodes 
at both poles because the corresponding basis function is 
Y 2 -1 oc k z (k x - iky). 

We also calculated the mass current for the 1/3- vortex 
and the line integral of the velocity field § vgdO along var- 
ious closed paths around the center. We find that the line 
integral yields approximately | when circling near the 
center, implying 1/3 quantization. However, this conclu- 
sion is not exact because the system is inhomogeneous, 
a situation different from the infinite system where the 
quantization is exact if the line integral is taken along a 
path far away from the vortex line. 
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